Hemoglobin signal network mapping reveals novel indicators for precision medicine

Precision medicine currently relies on a mix of deep phenotyping strategies to guide more individualized healthcare. Despite being widely available and information-rich, physiological time-series measures are often overlooked as a resource to extend insights gained from such measures. Here we have explored resting-state hemoglobin measures applied to intact whole breasts for two subject groups – women with confirmed breast cancer and control subjects – with the goal of achieving a more detailed assessment of the cancer phenotype from a non-invasive measure. Invoked is a novel ordinal partition network method applied to multivariate measures that generates a Markov chain, thereby providing access to quantitative descriptions of short-term dynamics in the form of several classes of adjacency matrices. Exploration of these and their associated co-dependent behaviors unexpectedly reveals features of structured dynamics, some of which are shown to exhibit enzyme-like behaviors and sensitivity to recognized molecular markers of disease. Thus, findings obtained strongly indicate that despite the use of a macroscale sensing method, features more typical of molecular-cellular processes can be identified. Discussed are factors unique to our approach that favor a deeper depiction of tissue phenotypes, its extension to other forms of physiological time-series measures, and its expected utility to advance goals of precision medicine.

Underlying the Precision Medicine Initiative is the premise that features of molecular expression are closely tied to disease subgroups 1 , and that applying individualized knowledge of these features will enhance understandings of disease susceptibility and treatment strategies.This focus has motivated the adoption of deep phenotyping strategies, of which multi-omics measures are considered a principal tool 2,3 .While informative, it has been recognized that often not leveraged by this initiative are varied forms of physiological recordings 4,5 .This is unfortunate, as this class of measures is recognized as information-rich as a consequence of sensitivity to a host of ongoing feedback adjustments by unseen drivers that support homeostasis and whose details are often disturbed by disease.While not directly observable, these unseen factors nevertheless produce macroscale behaviors affording the potential to identify surrogate signatures that might usefully serve to bridge macroscale behaviors with the unseen microscale phenomenology.
Invoked here is the rationale that identification of such factors can be gained through examination of lowvariance measures of short-term behaviors.An everyday example supporting this view is the seemingly hidden driving movements of a marionette when viewed from a distance.Composite adjustment of forces applied to varied strings (i.e., drivers) produce coordinated movements.Now consider observations that are time-delayed compared to the applied force.The greater the time delay, the more difficult it becomes to recognize the actions of drivers.It follows that reducing this delay will improve recognition of features causatively linked to such actions.
In a previous report we demonstrated that access to such low-variance measures can be gained by invoking a particular form of network analysis of time series 6 .Produced is a low-dimensional Markov chain whose recurring States are generated on a time scale similar to fleeting, nonstationary behaviors which, when combined with use of signal averaging over space-time, achieves the desired low-variance measures.When applied to hemoglobin time-series measures collected from women with breast cancer and a control group, revealed were multiple classes of adjacency matrices whose coefficient comparisons showed mainly uncorrelated patterns and disease sensitivity, suggesting that extended sensitivity to a host of unseen disease-sensitive drivers may be achievable 6,7 .
Here we have explored a logical follow-on to these measures and invoke the hypothesis that because multiple adjacency classes are generated, the actions of the considered drivers are likely dispersed across the same.Analogous to a puzzle, coalescing of coefficients across such matrices is taken as equivalent to assembling its

Transition-type patterns and disease sensitivities of network coefficients
Shown in Fig. 2 are selected group-mean adjacency matrices computed using the fixed-reference ordinal partition network (fr-OPN) approach described in Methods, applied to control-subject data.Plotted are results for State-based (Hb-State transition probability in Fig. 2a, pre-transition dwell times in Fig. 2b) and Hb componentbased measures (pre-transition mean values and fluxes for ΔtotalHb (Figs. 2c,e), and for ΔHbO 2 Sat (Fig. 2d,f)).Evident are varied contrasts compared to the absence of features in the cloud plots.Quantitative comparison of these features reveals correlations ranging from near zero (pre-transition dwell time vs. ΔHbO 2 Sat pre-transition mean) to a maximum of -0.68 (ΔHbO 2 Sat pre-transition mean vs. ΔHbO 2 Sat flux).The dissimilar patterns seen imply that varied information intrinsically exists in the measured time series but cannot be directly perceived from inspection of the primary data (Fig. 1b,c).Shown in Supplementary Figure 1 (Supplementary Note 2) is evidence that accompanying such varied patterns, are correspondingly dissimilar patterns of disease sensitivity.For the four network coefficients considered in Supplementary Figure 1, inter-breast comparisons yield p < 0.01 (unequal-variance t test) for 38-84% of transition types for breast-cancer subjects, but only 0-12% for noncancer subjects.
Whereas these results reveal different behaviors within the steady-state measures, possible systematic trends among two or more classes of matrices cannot be appreciated by inspection.Phenomenological reviews of such trends can be an efficient strategy for recognizing key drivers of system behaviors, which are often obscured in the case of complex systems.To accomplish this, we have systematically explored co-dependent behaviors by plotting the quantity in one adjacency matrix against that in another.In the case of paired comparisons, this produced a total of 153 plots, with additional measures involving three or more adjacency matrices also available.The binary plots explored all possible pairings among the 18 classes of measures: the three Hb-State (transition probability and pre-, post-transition dwell times), and 15 Hb-component (flux and pre-, post-transition mean amplitude for each of the five Hb-signal components) adjacency matrices.Examination of these co-dependency findings constitute the remainder of this report.

Evidence of structured co-dependent behaviors
If the transition-type patterns of the various adjacency matrices were as independent as inspection of Fig. 2 and Supplementary Figure 1 suggests, plots of one network coefficient against another would be expected to yield an unstructured distribution of points.While a minority of the plots generated do have this character (14%), most exhibit mainly smoothly varying trends including, unexpectedly, instances of fits to well-defined mathematical forms typical of saturable processes common in biology (e.g., enzyme catalyzed reactions).Because these findings represent a first report, where appropriate we have supplemented these observations with findings from modeling studies and noise measures, to document that trends seen have a biological origin and are not a consequence of applied mathematical transforms.
The ΔHbO 2 Sat versus ΔtotalHb plot in Fig. 3a (described more fully in Supplementary Note 3) is exemplary of a class of structured behaviors seen when the pre-transition mean amplitude for one Hb component is plotted against the pre-transition mean amplitude for another, for all 100 transition types.This figure reveals relationships that are not discernible by inspection of the same two coefficients when they are plotted separately (Fig. 2c,d).With 5 Hb-components, 10 paired plots can be generated, and the overall trends seen in Fig. 3a are obtained for all component pairings.To uniquely label each transition type, we applied different colors and symbol shapes and have further used open and filled symbols, as defined in Fig. 3b.The points in a given sector identify the paired pre-transition mean Hb-component amplitudes for the ten transition types that have a common pre-transition State.A similar plot of post-transition mean amplitudes for the same two Hb components would look almost the same as Fig. 3a, except that points in a given sector would have the same shape rather than the same color/ fill (Supplementary Figure 2, Supplementary Note 4).
It is instructive to consider whether aspects of the Fig. 3a predominantly linear patterning and observed sequence of points within a sector follow inevitably from the way that States are defined here and in Refs. 6,7.Our conclusion that they do not is established in Supplementary Note 5.The question was examined using the same model as adopted for the analysis of the State sequence, which ignores the radial dependence of data values (Fig. 1), in favor of an assumed uniform distribution.We computed the pre-transition mean-amplitude coordinates for all 100 transition types, under the null hypothesis that transitions between two points in ΔHbO 2 Sat-ΔtotalHb space will occur less often as the distance between the points increases.The findings were that the computed "spokes" are not linear, and that, unlike the pattern observed for real data, the post-transition State sequence differs across the pre-transition States, thus supporting the premise that the principal features of Fig. 3a are consequences of biological processes and not of any geometric dependence from applied State definitions.
The finding that trends in Fig. 3a depend on features of ΔtotalHb that are not otherwise evident (Supplementary Notes 3 and 5) raises the question whether this component or others influence trends seen in the individual adjacency matrices.To explore this, we implemented a strategy that compares parameter values for pairs of transitions that differ with respect to all Hb-signal components simultaneously, so that a dependence of .Represented on the color axis is the (occurrence percentage) 1/7 of temporally coincident resting-state ΔHbO 2 Sat and ΔtotalHb image values; the 7-th root, having a similar effect to logarithmic scaling, is plotted for ease of visualization.a given component-related parameter on any of them can be detected (Supplementary Notes 6 and 7).Applying this method, we find that, whereas the transition-dependent patterning of contrast features varies significantly between State-transition and Hb-component matrices (Fig. 2, Supplementary Figure 1), all have an underlying dependency on ΔtotalHb (Table 1, Supplementary Figure 3).While the magnitude of the bias varies, in every case where the sign of ΔtotalHb reverses (category '1' in Table 1) the average amplitude is reduced compared to when it is not changed (category '0').We further note that the magnitudes of these biases are greater than their counterparts for the other Hb-signal components (not shown) (p = 2.5 × 10 −10 -5.4 × 10 −3 (unequal-variance t test)).While physiological understandings may make the noted determinations unsurprising, they show that drivers of feature behavior can be effectively hidden.
Additional evidence of structured co-dependent behaviors can be gained from examination of plots of paired coefficient values before and after a transition has occurred.Shown in Fig. 3c is a plot of the same Hb components as in Fig. 3a, but here the post-transition mean value of ΔtotalHb is plotted on the x-axis instead of its pretransition mean.Inspection reveals that this substitution transforms the linear dependence to one having a welldefined hyperbolic functional form, a finding consistent with an underlying saturable phenomenon common among enzyme-driven behaviors.An alternative pairing shown in Fig. 3d involving substitution of ΔtotalHb flux on the y-axis, also reveals hyperbolic dependences, in this case extending to distinct trends for each pre-or post-transition State.As a qualitative measure of robustness of the data-value distributions in Fig. 3c,d, selected data points are annotated with two-dimensional (2D) error bars showing the standard errors for each plot variable.Unpaired and paired t tests were used for quantitative evaluations, as described in Supplementary Note 9.A.For the Fig. 3c pattern, the separations between data-point distributions for all nearest-neighbor pairs of pre-transition States are highly significant in one or both of the ΔHbO 2 Sat (p = 1.2 × 10 −9 -9.3 × 10 −4 ) or ΔtotalHb (p = 5.5 × 10 −5 -4.3 × 10 −4 ) dimensions.For Fig. 3d, all but the (largely superimposable) States 2,3 nearest-neighbor pair (and its reciprocal 7,8) have highly significant separations for both the pre-(p = 1.4 × 10 −11 -0.0097) and post-transition (p = 3.7 × 10 −11 -0.0029) groupings.
Examination of all 25 Hb-component pairings of flux versus post-transition mean value (Supplementary Figure 4) reveals that, while all produce smoothly varying trends, hyperbolic dependences are more limited than the linear trends seen for measures of mean component amplitude (Fig. 3a), being restricted to pairings that include the ΔtotalHb post-transition mean (Supplementary Note 8).Examination of the 25 pairings of pre-versus post-transition mean value likewise shows that hyperbolic dependences occur when ΔtotalHb mean values are plotted on one or both axes (Supplementary Note 8).For all such cases, an abbreviated representation is available from determination of the coordinates of the associated vertices and foci (Supplementary Figure 5, Supplementary Note 9) 11 .give the best-fitting approximations to the sets of same-shape symbols (i.e., fixed post-transition State).The curve-fitting function used 8 implements the generalized eigenvector fit algorithm of Ref. 9 .This method returns the coefficients for the best-fitting conic section: Ax 2 + Bxy + Cy 2 + Dx + Ey + F = 0.The conclusion that all fits were hyperbolic follows from the empirical fact that in every case considered the value of the discriminant (i.e., B 2 -4AC) was positive 10 .

Informative alternative groupings of transition types
More careful inspection of Fig. 3c,d reveals a pattern that suggests an alternative strategy for appreciating evidence of structured behavior.It is seen that all trends for a fixed pre-or post-transition State have a similar functional form, but each has a unique amplitude dependence (e.g., distance from the origin), thus identifying two distinct classes of behavior (i.e., a State-dependent amplitude within a State-independent functional form).Recognizing that a common dependence suggests sensitivity to a common underlying process, we explored an alternative classification of transition types into groups having a fixed dependence between pre-and post-transition States.www.nature.com/scientificreports/A straightforward way to define such groupings is to select all transition types that lie a fixed distance above or below the adjacency matrices' main diagonals (see inset of Fig. 4).These groupings are referred to as transition Classes.There are six Classes, numbered 0-5 in accordance with the distance, in State sectors, between the pre-and post-transition States.Each Class includes equal representations from all pre-and post-transition States, and all transition types in Class m (m = 0-5) have the property that exactly m Hb-signal components undergo changes of algebraic sign.
To demonstrate that transition Class information can provide additional insights into structured co-dependences, in Fig. 4 we have redrawn Fig. 3d with each marker relabeled according to its transition Class, revealing a systematic trend not evident by inspection of the original plot.It is seen that with increasing Class number, there are reciprocal trends in the ranges of values for the mean (decreasing) and flux amplitude (asymptotically approaches maximum value).The latter trend is present for all five Hb components (Supplementary Figure 6, Supplementary Note 10), demonstrating a binary response comprising a rapid rise phase (Classes 0-2) and an asymptotic phase (Classes 3-5).We conjectured that the bi-phasic trend reflects a transition from a larger or more heterogeneous set of unseen drivers affecting the Class 0-2 transition types to a smaller or more consistent set of influences on Classes 3-5.Accordingly, we anticipated that subsequent computations intended to reveal hidden-driver effects may yield notably different results for the 0-2 and 3-5 Class groupings.

Transition class-dependent enzyme-like behaviors in Hb-signal component fluxes
The findings of hyperbolic relationships among coefficients of the Hb signal components can be initially considered as reflective of unseen enzyme-like actions.To further explore this along lines of classical enzyme activity plots, we invoked analogies between transition probability and a hypothesized canonical substrate, and between (absolute value of) transition flux and an associated enzymatic velocity (Supplementary Note 11), and generated plots of 1/|ϕ k X | (Eq.( 10)) versus 1/P k (Eq.( 5)) for each Hb-signal component.Inspection reveals suggestions of (curvi) linear structure, but also notable deviations from the structured behavior.This is illustrated in Fig. 5a, for the case of X = ΔtotalHb and ϕ k X computed using the SMTM averaging approach (Supplementary Note 12).A classical Lineweaver-Burk (L-B) analysis was performed that was limited to Classes 3-5 data values (black symbols in Fig. 5a inset).While the result is statistically significant (r = 0.38 (p = 0.018)), the dependence of 1/|ϕ k X | on 1/P k is complex.Nevertheless, a large number of Classes 3-5 points plotted do show dependence consistent with a saturable process, with the deviations limited to a high transition-probability subset and nearly all elements of transition Classes 0-2.A second L-B computation, which considers only the non-deviating points (n = 38), yields a highly linear dependence (i.e., dashed curve in Fig. 5a inset; note that ordinate data values are log transformed), with r = 0.989 (p < 10 −10 ).We tentatively concluded that an enzyme-like relationship is present, but it is only one of the factors that determine the overall of 1/|ϕ k X | versus 1/P k dependence.Examination of the full range of Hb components, averaging methods, and breast groups revealed that the subset of transition types that participate in the structured dependence is different, and seemingly arbitrary, for every component.This is illustrated in Fig. 5b,c, where Fig. 5b is a plot of 1/|ϕ k X | versus 1/P k for X = ΔHbO 2 Sat (r = 0.90 (p < 10 −10 )), while in Fig. 5c the dependent variable is the reciprocal of (the absolute value of) the "dwelltime flux" defined as the difference between post-and pre-transition dwell-time values (Supplementary Note 13) (r = 0.94 (p < 10 −10 )).The composite of findings was taken as unsurprising, since flux values are expectedly influenced by multiple factors, but also as encouraging, in that the (assumed) enzymatic activity factor is sufficiently strong as to reveal itself in the plots at all.
Based on the parameter-sweep findings, we modified the double-reciprocal plotting procedure in a way that yielded improved evidence of structured behavior: flux values for each subject were transformed to a normalized quantity as described in Methods (Eq.( 11)), which reduced the impact that inter-subject variability in flux amplitudes has on plotted group-mean values.Consequently, linear scales can be used for both axes, and there is a notable disease dependence.However, the Class 3-5 data points are still not aligned linearly and no nonarbitrary procedure is discernible for excluding the non-aligned transition types for any result based on a single Hb component.
The insight that enabled further development of the L-B analysis came from the recognition that the specific transition types that deviate from the structured behavior are different for each component considered in Fig. 5a,b,c.This raises the possibility that composite measures obtained by vectorially combining flux data for multiple components (in particular, the t′-scores, which are dimensionless, thereby avoiding incompatibleunits concerns) may reveal a structured behavior that includes all Class 3-5 transition types.Shown in Fig. 5d is the result obtained when the mentioned refinements were applied to a 3-component flux measure (Eq.( 11)) (Supplementary Note 13).Here we find an excellent linear fit over the entire span of coefficient values, demonstrating that off-trend values in each one-component fit are overcome by contributions from the other components.There are highly significant differences between both the Fig. 5d regression slopes (p = 1.6 × 10 −5 ) and the y-intercepts (p < 10 −10 ), and selected Fig. 5d data points are annotated with 2D error bars (± 1 standard error for each plot variable), as a qualitative measure of robustness of the linear trends.Similar-quality fits were obtained for a selected individual (slopes: p = 0.0026; intercepts: p = 1.1 × 10 −6 ; Supplementary Figure 7), emphasizing that neither the qualitative aspects nor the goodness of fit require variance reduction from group averaging.Notably, the effect in tumor subjects is to enhance V max , which is consistent with the known action of NO, whose levels are increased in breast cancer 12 , on the enzyme soluble guanylyl cyclase 13 , a key driver of vascular smooth muscle contractility.We invoke this interpretation in recognition that vascular pulsatility is the dominant phenomenology affecting Hb signal dynamics.www.nature.com/scientificreports/

Evidence for tumor biomarker sensitivity in L-B plot coefficients
For the Fig. 5d demonstration of disease sensitivity, the plotted T-group vector amplitudes are averages over all breast-cancer subjects.We have extended the analysis to investigate whether the vector-computation results exhibit sensitivity to either the ER or Her2 breast-cancer biomarkers, by performing separate L-B analyses for the ER( +) and ER( −) subjects (n = 13 and 5, respectively), and for the Her2( +) and Her2( −) subjects (n = 9 and 9).Statistically significant differences between the regression slopes were found in both cases (ER: p = 2.2 × 10 −8 , Her2: p < 10 −10 ), and between the ER( +) and ER( −) y-intercepts (p = 0.0081) (Supplementary Figure 8, Supplementary Note 14).Additionally, t tests were applied to vector amplitudes averaged across subjects (i.e., one net value per transition type) to compare the amplitudes for marker-( +) and marker-( −) subjects.Tests performed considered either all Classes 3-5 transition types, or only those for one of the four vector-component rankings identified in Table 2. Results demonstrate that vector amplitudes are sensitive to the presence or absence of a

Discussion
The understanding that the varied organizational scales of biology and their homeostatic balance are strongly influenced by lifestyle, environmental factors and varied forms of molecular expression has motivated a focused buildout of resources to achieve a granular description of molecular expression and their responses to principal drivers 1,14,15 .This effort builds on the idea that descriptions of this type will support similarly individualized approaches to disease prevention, early detection, stratification and therapeutic interventions.Recognized as precision medicine, its implementation thus far is mainly built on gaining deep phenotypic descriptions of the compositional elements of biology.Generally, not considered by this information mix are varied forms of physiological measures 4,5 .Reflecting the integration of a host of body functions, the relationship between macroscale behaviors recognized by such measures and molecular-scale descriptions and associated phenomenology are restricted principally to qualitative trends.As such, the information densities achievable from physiological measures are orders of magnitude lower than multi-omics measures.
It is our view that this dichotomy is not due to any fundamental information barrier, but rather speaks to current limits on feature recognition.The understanding that drivers of time-varying phenomenology are hidden to the sensing method motivates adoption of methods that can serve as useful surrogates.Among the varied approaches used to explore dynamic behaviors, attention is given to network methods with the understanding that useful descriptions of system behaviors can be appreciated without requiring detailed knowledge of the drivers of system dynamics.Experience shows that network formulations of physiological time-series measures (e.g., hemodynamic, bioelectric) is typically directed to assessing a network's structure and stability 16 .While features of system organization are thereby revealed, network measures that explore system dynamics per se have been mainly absent.As recognized, the common approach of reducing information to a single metric (e.g., temporal correlation) strongly limits the granularity of accessible information compared to the rich details of the native signal 17 .Lost are dependences between the fleeting originating signals, whose short-term amplitudes vary owing to the feedback actions of a multitude of unseen factors.
Being data-driven, State descriptions are free to consider any of a host of observable phenomena.Here we have extended principal features of the hemoglobin signal (i.e., oxyHb, deoxyHb) to include other computable quantities (i.e., HbO 2 Sat, totalHb, and HbO 2 Exc).As the latter are dependent phenomena, it might be thought that they are simply redundant information.However, their inclusion provides a principled approach to objectively expanding the definable network space (Supplementary Note 1).Also, as recognized, the State definitions applied here are coarser than necessary 6 .Descriptions that are more restricted in the spatial and time domains are available, including the history dependence of State transitions.In the limit, network descriptions can be generated from even a single measured time series.Also available are descriptions that make use of information acquired from multiple forms of simultaneous measures.In terms of the hemoglobin signal, attention is drawn to the technique of diffuse correlation spectroscopy 18 , which provides measures of blood flow and has been shown sensitive to elements of the tumor phenotype 19 .Inclusion of these would extend the network to six components, rather than the five explored here.Similar considerations would also apply to other forms of timevarying physiological measures (see below).
Table 2. Impact of tumor marker status on tumor-breast flux-vector amplitude.Italics indicate statistical significance (p < 0.05).The tabulated rank-order categories refer to the relative magnitudes of the three normalized quantities (i.e., t′-scores) used in computing vector amplitudes (Eq.( 11)).Specifically  www.nature.com/scientificreports/Further recognized is the particular approach used to define network features.Here we assume that signals satisfying the same State definition are responses to similar levels of underlying drivers.It is our contention that this assumption is more strongly justified under the fixed-reference network-node definitions used in this report and in Ref. 6 (fr-OPN; see "Methods") than would be true if the time-varying multivariate OPN method 20 had been employed instead (Supplementary Note 17).As outlined, the former supports more granular descriptions of network behaviors.Such noted flexibility of State descriptions, combined with other sensing measures, suggests that access to a broad range of small-scale behaviors might prove feasible.
As alternative network to time series methods are available 17 , it is useful to consider their expected utility for the aims of this study.Recognized are three principal strategies comprising use of either recurrence, visibility, or ordinal partition operations (RN, VN, and OPN respectively).Comparisons among them directs attention to how network nodes and edges are defined.For univariate measures, RN or VN methods are less desirable because network parameters are defined from sequences of temporal information, thus convolving temporal behaviors.Similarly, for RN, while multivariate methods can limit temporal convolution in node descriptions, a comparable impact on the variable recurrence times associated with network edges is lacking.Also, as mentioned, while more typical implementations of multivariate OPN methods can be used, these have the effect of blurring assignments of network transitions relative to those considered here (Supplementary Note 17).Overall, because the goal of these methods has been to appreciate topological features, it is unsurprising that their typical uses are not well suited for the aims of this study.It follows that appropriate exploration of other phenomenology (e.g., short-term dynamics, thermodynamics) can be expected to require careful curation of network properties.
It is likewise appreciated that well-established non-network time-series methods, while useful for classification 21 or feature-extraction [22][23][24] , and potentially useful for pre-conditioning operations for the method of this report, do not generate data from which low-variance measures of short-term dynamics can be directly read.Additionally, while these methods lie on a spectrum with regard to the requisite degree of prior knowledge of the processes governing system dynamics (e.g., low for wavelets 21,22 , high for Kalman filtering and related methods 25,26 ), mathematically they are more strongly assumption-dependent or model-based than the fr-OPN method used here.
Extension of the applied methodology to other forms of physiological monitoring also appears feasible.For the case of EEG or MEG, time-varying amplitudes of different frequency bands can serve as a mathematical counterpart to the hemodynamic components.The time series for each band would have a temporal mean value, with the instantaneous amplitude alternately above ( +) and below ( −) it, allowing for the definition of States and adjacency matrices in a manner directly parallel to that used for the fNIRS Hb-signal components.While adjacency-matrix computation could include temporal integration over the full measurement period as in this report, also plausible would be combining the time series of frequency-band States with those for EEG microstates computed as described in Ref. 27 .For example, the sub-intervals corresponding to only a specified microstate could be selected, and these contiguously arranged to create a modified time series.The Statenetwork analysis applied to microstate-specific time series could help distinguish between intra-microstate and microstate-switching phenomenology.
For fMRI, counterparts to the EEG microstates can be computed from resting-state image time series 28 , but equivalency to the EEG frequency bands is mainly lacking due to lower sampling rates.As an alternative, a States definition can be based on the algorithmic feature that, for every image frame, a probability of each microstate being the "true" one is generated 28 .For their purposes, the authors took the reasonable approach of changing the largest probability value to 1 and all others to 0. To define counterparts of the Hb-signal components and States, we would propose instead to compute mean values for each microstate's probability time series, and at every time frame note which microstates have probabilities above ( +) or below ( −) their temporal means.In this way, each microstate (which is a spatial distribution of BOLD-signal amplitudes) is regarded as a signal component, and the States as different combinations of probability-based algebraic signs.From the State and componentamplitude time series, network adjacency matrices would be computed in the same manner as reported here.
Returning to the goal of including physiological measures as indicators for precision medicine, our demonstration of specific coefficient patterning across the varied adjacency matrices (e.g., Fig. 2 and Ref. 6 ), sensitivities to ER and Her2 expression, and the ability to recognize particular forms of structured behaviors (e.g., Fig. 5d), supports the hypothesis that the considered methodological extensions would likely identify associations to other forms of small-scale behaviors.In one form, recognition of these could look to parallel methods currently applied to explore multi-omics findings to identify correlations between feature classes and their dependence on varied disease characteristics 15 .For fNIRS measures, attention is drawn to inflammatory responses and factors impacting oxidative metabolism.Another would be to recognize dependencies between specific pharmacological manipulations and their impact on individual coefficient patterns or their co-dependencies.Then, because fNIRSbased network features and multi-omics information reflect different biological determinants and time scales, in combination it may be possible, for example, to improve on the reported accuracy of the PINSPlus-tumor subtyping algorithm 15,29 .While we have emphasized recognition of structured behaviors, more generally it seems likely that disease-sensitive coefficient patterns need not follow any particular mathematical form.In such cases, exploration of coefficient matrices using varied data classification schemes would be warranted.
In summary, we have shown that by invoking methodology that can reveal low-variance measures of shortterm dynamics, and by exploring their higher-order co-dependencies, evidence of behaviors consistent with actions arising from molecular-cellular processes is revealed, even though the primary sensing method is macroscopic.We expect that this methodology can be extended to other forms of physiological monitoring.We further recognize that extension of State class descriptions and their expected history dependences will aid in identifying the actions of other forms of unseen drivers which, when combined with multi-omics methods, will extend the information base upon which decisions of precision medicine depend, while offering the advantage of on-demand monitoring.www.nature.com/scientificreports/Summing U k over its spatial dimension and transposing yields C k , the 1 × (N t -1) matrix of type-k transition counts, whose i-th element is: Performing the computation in Eq. ( 3) for all 100 transition types produces C, a 100 × (N t -1) matrix of transition counts, where C k is the k-th row of C.
As an illustration of method, we show an example of a C matrix in Supplementary Figure 9.While the qualitative properties of the plotted C matrix (e.g., alternately high and low counts), are typical, the durations and starting times of the transition-count bands are subject-specific.Thus the information in C is not immediately amenable to inter-subject comparison and quantitative analysis.We convert it to a more tractable form by, first, summing over the temporal dimension to produce the 100-element transition-count vector c, whose k-th element is: and N TS = N t -1 is the number of measurement time steps.Each c k value is directly proportional to the measurement duration, which ranges from 400 to 1000 time frames for different study participants (median value and mode = 700, and disease and non-disease subject groups are measurement duration-matched (p = 0.14)).Accordingly, for inter-subject comparisons we normalize c to the sum of all values in c: The k-th element of (dimensionless) P is the probability for the k-th transition type, in units of percent.
It should be noted that the summation order in the preceding description-first across image voxels (Eq.( 3)), then time (Eq.( 4))-could be reversed, or the U k matrix (Eq.( 2)) could be summed across both dimensions simultaneously, without impact on the computed values of transition probability (Eq.( 5)).In contrast, computed values of the network parameters defined in the remainder of this section are dependent on the sequence of operations.Simultaneous summation across both dimensions-the "grand average" (GA)-is the method that is most straightforward, is presented in subsequent equations, and was used to generate all Results in this report except where indicated otherwise.A description of the alternative data-averaging schemes-"temporal mean of the spatial mean" (TMSM), "spatial mean of the temporal mean (SMTM)-and a consideration of circumstances in which their use may provide useful counterpoints to GA, is presented in Supplementary Note 12.

Computation of inter-transition dwell times
There are ten transition types that correspond to a voxel being in the same Hb State in successive time frames (using the single-number indexing defined in Eq. (1), they are types 11n -10, where n = 1-10).The term adopted here to describe occurrences of these transition types is that the voxel "dwells in" State n during that time step.For each transition type k, two distinct average dwell times are extracted from transition matrix T. These are the mean number of time frames that voxels dwell in the pre-transition state ( s 1 = ⌈k/10⌉ , where ⌈x⌉ (i.e., the ceiling function) is the smallest integer ≥ x) prior to transition k, and the mean number of frames that they dwell in the post-transition state (s 2 = k -10(s 1 -1)) following it.We introduce the quantities n k s 1 ij and n s 2 k ij , to denote the pre-and post-transition dwell times, respectively, for the i-th time-step transition in the j-th voxel.That is, n k s 1 is the number of time frames dwelt in State s 1 prior to a type-k (i.e., from s 1 to s 2 ) transition, and n s 2 k is the number of time frames dwelt in State s 2 following a type-k transition.We further introduce τ k (1) and τ k (2) to denote the mean pre-and post-transition dwell times, respectively, for type-k transitions.Thus, we have: where (U k ) ij (Eq.( 2)) is 1(0) if a type-k transition does (does not) take place in the j-th voxel in the i-th time step, and c k (Eq.( 4)) is the total number of type-k transitions in the image time series.The two dwell times have different time-step summation limits in the index-i summations of Eq. ( 6), because pre (post)-transition dwell time cannot be evaluated for transitions that occur during the first (last) time step.
For the ten s n → s n transition types there can be two or more consecutive occurrences of the same type, and this would produce multiple values of n k s 1 and n s 2 k for a single dwell period.In these instances, only the largest values of n k s 1 and n s 2 k in each sequence of consecutive occurrences are included in the Eq. ( 6) summations.

Quantification of transition-linked hemoglobin-concentration and -saturation changes
The Hb States and transitions are defined in a manner that considers a fixed-reference ordinal dependence (i.e., independent of component amplitudes, and hence of distance from the axes or origin in Fig. 1).At the same time, image time series contain information on the magnitude of each Hb-signal component, which can be combined www.nature.com/scientificreports/with the categorical information to compute the average levels of the five components of the Hb signal before and after each type of transition, as well as the average amounts by which those levels change during transitions.The starting point for the concentration-and saturation-dependent quantities is five matrices, one for each component of the Hb signal, defined in a similar manner to the T matrix in Eq. (1): where ΔD, ΔE, ΔO, ΔS and ΔTo denote the image time series (formatted as N t × N v matrices) of ΔdeoxyHb, ΔHbO 2 Exc, ΔoxyHb, ΔHbO 2 Sat and ΔtotalHb, respectively, while Δ 2 D, Δ 2 E, Δ 2 O, Δ 2 S and Δ 2 To are the corresponding concentration and saturation changes that attend each transition.
Proceeding in parallel with the transition probability derivation, the analogue of Eq. ( 2) is: Thus, the arrays with superscripts '(1)' , '(2)' and '(3)' contain pre-transition, post-transition, and transition-associated change data for the k-th transition type, respectively.We sum the arrays defined in Eq. ( 8) over position and time, and account for inter-subject variations in measurement duration by normalizing to the transition count: and The (μ k X ) (1) parameter is the mean value, per transition, of pre-transition concentration or saturation, while (μ k X ) (2) is the corresponding post-transition value.We refer to ϕ k X as the flux for signal component X and transition type k. (A related quantity, called the transition mass in Ref. 6 , is computed by substituting the product of N TS and N v for c k in Eq. (10) but is not considered in this report.)

Computation of higher-order adjacency matrix co-dependent values
If a selected number M of State-or transition-derived parameters are treated as coordinate axes of a mathematical space of M dimensions, then the parameter values for the k-th transition type can be plotted as a single M-dimensional point in that space.For most of the results presented in this report, M = 2 and the parameter values are plotted in their original, dimensional units (e.g., seconds for dwell times, molar (concentration) for ΔtotalHb pre-transition mean values).However, also explored is a variation on the same idea, in which the k-thtransition-type parameter values are converted to dimensionless units; this additional, normalization step allows us to compute an amplitude (and direction, if desired) for the resulting M-dimensional vector.
A particular example considered in this report is the vectors obtained for the three-parameter combination of ΔtotalHb transition flux, ΔHbO 2 Sat transition flux, and transition dwell time.A vector amplitude cannot be directly computed as A = [(ΔtotalHb flux) 2 + (ΔHbO 2 Sat flux) 2 + (Δ(dwell time)) 2 ] 1/2 (where Δ(dwell time) = (posttransition dwell time) -(pre-transition dwell time)), because the quantities considered have different units.To eliminate unit inconsistency, each parameter was converted to dimensionless t′-scores: t k ′ = (x k -m)/s, where x k = parameter value for the k-th transition type, m = mean value across all 100 transition types, s = standard deviation across all types.However, in order to preserve inter-breast disparities that may be present in the original data 6 , the data values for each breast were referenced to the mean and standard deviation for the contralateral breast (e.g., t′(T) = (x(T) -m(U))/s(U), where T, U = tumor-bearing and contralateral unaffected breast, respectively, of subjects who have unilateral breast cancer; t′(L) = (x(L) -m(R))/s(R), where L, R = left and right unaffected breast, respectively, of non-cancer subjects).Vector amplitudes for the k-th transition type and breast group X (where X is one of L, R, T or U) were subsequently computed as:

Study participants
Findings reported herein were derived from publicly accessible, de-identified data from 63 subjects-18 with confirmed breast cancer (6 right-breast, 12 left-breast, average tumor size ~ 2.7 cm, range 0.5-6 cm), and 45 age-and BMI-matched control subjects, 23 of whom had evidence of non-malignant pathologies in one or both breasts 35 .A more detailed description of enrolled subjects is given in Ref. 30 and in Supplementary Table S1.

Figure 1 .
Figure 1.Mapping of continuous time-series signals into discrete Hb states.(a) Each Hb State is a region corresponding to one of the ten physiologically allowed permutations of algebraic signs of the five components of the Hb signal, depicted as regions in the ΔtotalHb-ΔHbO 2 Sat coordinate system (circled numbers are State assignments).Each secondary axis is the locus of points at which one of the other three components is equal to its baseline mean value (e.g., ΔdeoxyHb = 0 on the positively sloped axis).Colored arrows indicate which regions of the plane correspond to positive (solid-line arrows) and to negative (dashed-line arrows) deviations.(b-c) Paired resting-state ΔtotalHb and ΔHbO 2 Sat values for affected (b) and unaffected (c) breasts of a selected breast-cancer subject (50 y/o, BMI = 44, 4 cm intraductal carcinoma).Represented on the color axis is the (occurrence percentage) 1/7 of temporally coincident resting-state ΔHbO 2 Sat and ΔtotalHb image values; the 7-th root, having a similar effect to logarithmic scaling, is plotted for ease of visualization.

N 2 .Figure 4 .
Figure 4. Transition-Class Trends Embedded in Binary Co-Dependences.Replot of the ΔtotalHb flux versus mean post-transition ΔtotalHb data in Fig. 3d, with points labeled according to the transition Class corresponding to each transition type (inset).Dashed line segments show the maximum values of ΔtotalHb flux (horizontal segments) and mean post-transition ΔtotalHb (oblique segments) for each Class, demonstrating an inverse relationship between the flux and mean-value coefficients. https://doi.org/10.1038/s41598-023-43694-7

Table 1 .
Distribution of adjacency-matrix values between transitions with/out reversal of ΔtotalHb algebraic sign.Italics indicate statistical significance (p < 0.05).a '1' ('0') denotes transition types for which the algebraic sign of ΔtotalHb does (does not) change.
bUnequal-variance unpaired t-test; n('0') = 40 (adjacency-matrix main diagonal is excluded), n('1') = 50.c Breast group N is the union of the U, L, and R groups.d Correlatedsamples t-test; inter-group correlations range from 0.82 to 0.998.e Input data are averages of the pre-and posttransition means.